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Abstract. We present an exposition of a method of discretizing ordinary differential 
equations while preserving their Lie point symmetries. This method is very general and can 
be applied to any ODE with a nontrivial symmetry group. The method is applied to obtain 
numerical slutions of second and third order ODEs invariant under two different realizations of 
S L(2, R). The symmetry preserving method is shown to provide a better qualitative description 
of solutions than standard methods. In particular it provides solutions that are valid close to 
singularities and beyond them. 
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1. Introduction 

Lie group theory was originally invented as a systematic tool for obtaining exact analytical 
solutions of ordinary and partial differential equations(ODEs and PDEs). For ODEs the ex- 
istence of a nontrivial symmetry group (local Lie group of local point transformations taking 
solutions into solutions) makes it possible to reduce the order of the equation. If the symme- 
try group is large enough, the problem can be reduced to quadratures [|T1. The reduction to 
quadratures in principle provides the general solution of the ODE. This may however be in 
implicit form that is of little use in visualizing the solution, or presenting it in the form of a 
graph. What happens in such cases is basically that an ODE is transformed into an algebraic, 
or transcendental equation. 

If the symmetry group is not large enough, or its structure is such that the solution pro- 
vided by the group is too implicit to be useful, it is necessary to resort to numerical solutions. 
The question arises: can the symmetry group still be put to good use? 

All numerical methods for solving ODEs replace the differential equation by a difference 
one, usually on an a priori chosen lattice, either a regular one, or one adapted to some known 
or expected behaviour of solutions. Most or all Lie point symmetries are lost in this procedure. 

Over the last 20 years a considerable effort has been made to apply Lie group theory 
to difference equations (for reviews containing references to the original articles see [[21 0). 
One approach to this problem is to generate invariant difference schemes, consisting of two 
equations, involving independent and dependent variables (jc„, j„) evaluated at A'^ -I- 1 points 
(in order to be able to approximate an ODE of order A^^) 

Ea{Xn+K,-,Xn+L,yn+K,-,yn+L) =0, a = 1,2, L-K = N -\. (1) 

The equations must allow to calculate xn+l, Jn+l if all previous points are known. This 
is a condition on the Jacobian matrix, namely 

.aJ- -) = 2. (2) 

Thus the scheme ^ represents a difference equation and an equation for the lattice. 

Let us consider a Lie algebra g realized by vector fields in two variables 

Xi = ^i{x,y)dj, + (pi{x,y)dy i= l,...,M = dim{g) (3) 

where x and y are respectively the independent and dependent variables in a differential, or 
difference equation. 

Invariant ODEs of order A^^ are obtained as follows. We prolong the vector field (|3]) to 
order N 

pr^Xi = Xi + (P^dy + (P^'dy' + ... + (f>f'dyiN, , (4) 
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where the coefficients 0^ , are expressed in terms of partial derivatives of the coefficients 
^, in the original vector field (|3]) flU. Invariants /^(x, —,yNx), 1 < H < A, of the group 
action corresponding to Lie algebra g are obtained by solving the system of first order PDEs 

prXi(p{x,y,y^, ...,yNx) = 0, / = 1, ...,M. (5) 

The number of functionally independent solutions of (|5]) is at least A = A'^ + 2 - M. It can be 
larger if the M equations ([5]) are linearly dependent on some manifold S in the corresponding 
jet space. The invariant equation will have the form 

F(/i,...,/^) = (6) 

where F is any sufficiently smooth function. 

An invariant difference scheme is obtained in a similar manner. We write the vector fields 
(|3]) in one chosen point n then prolong them to as many points as figure in the scheme (i.e. 
+ 1). We have 

L 

pro Xi^n = ^ [^i(Xn+k,yn+k)dx„^, + (pi(Xn+k,yn+k)dy„^t] I < i < M. (7) 

k=K 

The discrete invariants In,a(Xn+j,yn+j), 1 < a < B, K < j < L axe obtained by solving the 
system of equations 

pro Xi^n(f){Xn+j, yn+j) = 0. (8) 

The number of functionally independent solutions of ([8]) is at least B = IN + 2 - M and is 
larger if the equations dsl) are linearly dependent on some manifold § in the discrete jet space. 



In the continuous limit (|7]) reduces to (|4]) so the discrete invariants will reduce to the 
continuous ones. In general there are more discrete invariants than continous ones, so some 
combinations of the discrete invariants will go to zero in the limit, others will go to the 
continuous invariants. This makes it possible to choose an appropriate basis for the discrete 
invariants and to write an invariant difl'erence scheme: 

£i =F(/„,i,...,/„,a) = (9) 

E2 = E2(In,U...,InA) = (10) 

with F as in (|6]) and E2 satisfying £2 ^ in the continuous limit. 

Such diff"erence systems may be of interest in their own right and describe discrete phe- 
nomena on some specific symmetry adapted lattice. On the other hand the difference scheme 
may be chosen to have a specific ODE as its continuous limit. By construction the ODE 
and the difference scheme will be invariant under the same symmetry group G. Solving the 
difference scheme numerically provides approximate numerical solutions of the ODE. Since 
the symmetry group G determines many properties of the solution space one can expect that 
numerical schemes using a symmetry adapted discretization will have some advantages over 
other numerical methods. It has indeed been shown that for first order ODEs symmetry pre- 
serving discretizations are exact: the invariant differential equations and difference schemes 
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have exactly the same solutions H. Symmetry preserving discretizations of second order 
ODEs can be solved exactly using a Lagrangian approach (51161. These analytic solutions of 
the difference schemes then converge rapidly to the solutions of the ODEs [6J. Two recent 
articles dTk iSJ were devoted to numerical solutions of second and third order ODEs. It was 
shown (at least for the considered examples) that the qualitative behavior of solutions of the 
ODEs, specially in the neighbourhood of singularities, is better described by symmetry pre- 
serving discretizations than by standard methods. 

Four inequivalent realizations of sl{2,R) by vector fields of the form ([3]) exist BH. In 
this paper we concentrate on two of them, not treated in previous articles [|7lj8l. We construct 
their diff"erential invariants up to order three and their diff"erence invariants involving up to four 
points. This allows us to write all invariant ODEs of order up to three and their discretizations. 

In Section 2 we present the four realizations of sl(2,W). The invariant ODEs and their 
discretizations are presented in Section 3. Section 4 is devoted to numerical solutions and 
Section 5 to conclusions. 



2. The four realizations of sl(2,R) 

Let {Xi,X2,Xi,} be three vector fields of the form (|3]) satisfying the commutation relations 

[XuX2]=Xu [X2,X,]=X,, [XuX,]=2X2. (11) 

One of them, say Xi can be straightened out to Xi = dy. Then X2 can be transformed either 
into X2 = ydy or X2 = xd^ + ydy. Point transformations leaving the standardized fields Xi and 
X2 invariant will further simplify X3 and we obtain the four inequivalent realizations, namely: 



1. ^/i(2,R): 

X,=dy, X2=ydy, X3=y%. (12) 



The three vector fields (12) are linearly connected, that is in any given point of they are 
linearly dependent. This sl(2, R) algebra is not maximal among finite dimensional subalgebras 
of diffi2, R) but can be imbedded into slj,(2, R) e sly(2, R) with 

sl,i2,R) = {d_„xd„x'^d,,}. (13) 

For the remaining three 5/(2, R) algebras no two of the three vector fields Xi, X2 and X3 are 
linearly connected. 



2. 5/2(2, R): 

X, = dy, X2 = xd, + ydy, X, = 2xyd, + y^dy. (14) 

This 5/(2, R) algebra is not maximal in diff{2,W). We can add X4 = xdx and obtain the al- 
gebra gl(2,W). The algebra is imprimitive in that the coefficients of dy are aU functions of y 
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3. 5/3(2, R): 

X,=dy, X2 = xd,+ydy, X, = 2xyd, + (-x'+y^)dy. (15) 

4. 5/4(2, R): 

Xi = dy, X2 = xd, + ydy, X3 = 2xyd, + (/ + (16) 

These two realizations are equivalent over C but not over R. They are both primitive and both 
are maximal subalgebras of diff{2, R). 

The ODEs invariant under 5Li(2, R) and 5'L2(2,R) were treated earlier [|7l[8l. Here we 
concentrate on SLj,{2, R) and 5L4(2, R). 

3. Invariant ODEs and difference schemes 

3.1. 5/3(2, R).- Xi = dy, X2 = xd^+ydy, X3 = 2xyd,, + (/ - x^)dy 

A complete set of functionally independent differential invariants up to third order is 

^ ya+y^)-xr ^ 3xyy"'-xY\l+y'') 

1 (l+y2)3/2 ' 2 (i+y2)3 • ^ ^ 

Note that a complete family up to any order can then be deduced using P. Olver's Propo- 
sition 2.53 stated in Section 2.5 of llj. This is true for all realizations. 

In the discrete case a basis for all 3 point invariants is 

jn _ ( (^n - Xn-lf + (Vn - yn-\f ^'^ _ I (Xn+l - Xnf + (jn+l - ynf ^'^ 

\ XfiXn-\ I \ Xfi^iXfi I 

jn+l _ ( (^n+i - ^n~\f + (jn+l ~ yn-l)^ \ ^ ^^^^ 
\ Xn+ 1 Xfi- 1 / 

A complete family for any number of points can then be obtained simply by shifting the 
above invariants. Namely, the shifts of /"^^ and of I'^'^^ would be the two new invariants for a 
5 points scheme. 



Combinations of the discrete invariants ( 18 1 that approximate /i and I2 from ( 17 1 are 



'1 ' "1 

respectively. 



™+i _ _g J ^1 1 ^ + 1 and (19) 

^ \ H^r +^1) / 

= — 1 — T f^r" - jv') 
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(i+y2)3/2 (20) 

where C is an arbitrary constant. To get rid of possible sign ambiguities, we solve the square 
of ( 20 1 /j = to obtain the solution 

{y-yof + ix±C/af = l/a^ (21) 

with a,yo integration constants and a 0. Those are circles with center (+C/a, jo) and radius 
r = 1 1 a. 



An OaS that goes to the ODE (20) in the continous limit is obtained if we put 

jn^i = _8 J — /I 1 +1 = c, E2{riir\rr^) = o iii) 
^ \ i"j'i+\ri+^ + II) I iv 1' 1 ' 2 ^ 

with E2 defining the mesh and going to in the continuous limit. 
The 3rd order invariant ODE can be written as 

3xVj"' - ^V"(i + /') _ ^ l y'a+y'^)-xy" \ ^23) 



(l+y2)3 ^ (l+3;'2)3/2 



An OaS that goes to the ODE (23 ) in the continous limit is obtained if we put 



= T(-/r' - J^^') = F{JT') (24) 

-r -r ij 



where 7"^ is given in (22 1 and the lattice is 



E2ii"„ir\ir',i2"\i'r') = o (25) 

with E2 going to in the continuous limit. 

3.2. sk{2, R): Xi = dy, X2 = xd, + ydy, X3 = Ixyd^ + {x^ + y'^)dy 

A complete set of functionally independent differential invariants up to third order is 
^ xy"+/(jy'2-l) 



h = 



(y2 _ 1)3/2 ' 

Ix^iy' + 1)/" + 3((/ - \){y' + 1)^(3/2 - 1) + Axy'iy' + 1)/' - 2x^y"^) 



(26) 



(y' - l)2(y + 1)3 
and the second and third order ODEs are 

h = C, (27) 
h = F(h). (28) 
We again take the square of the second order equation = and obtain the solution 

(x±C/af-(y-yof = l/a^ (29) 
with a 0,yo integration constants. The solutions for a are hyperbolas. 
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In the discrete case a complete set on 3 points is given by 

^2 V- ^2 \i/2 



rn+l 

' I - 



rn+l _ 
^2 - 



4XnX„^l - iiy„ - J„_i)2 - (Xn - X„„i)2) 

AXn+\Xn-~\ - {iyn+\ - Jn-lY - {^n+l - Xn-\Y) 



(30) 



1/2 



Combinations of the discrete invariants (30) that approximate I\ and h from (26) are 

1/2 



,..1, V^( /2-(/i+/i.) _^ 



(31) 



= — r(7r" - rr) + 67f + 3 

-r -r ij 

SO the corresponding OaS are respectively 
J"^-"^ = C, 7f 2 = ^(7;^+') 



(32) 



where C and F are the same as in (27 ) and (28 ). An invariant equation for the mesh must be 



added in each case as in (22) and (25 ). 



4. Numerical solutions 

To perform numerical tests, the arbitrary function F of Section 3 needs to be specified. We 
choose F(/i) = (this choice is arbitrary). 

The second and third order ODEs are then discretized in several different ways: standard 
finite diff"erence methods, matlab solver (ode45) which uses Runge-Kutta methods of order 
four and five, and the symmetry preserving discretization which has been described in the 
previous sections. The standard finite difference method consists of approximating the 
derivatives of the dependant variable using Lagrange interpolation polynomials (see ^ for a 
more detailed explanation of the numerical methods used). As a quick reminder, the standard 
finite difference method gives on a four point scheme 

^ -(27(j„+l - yn) - {yn+2 - yn-\)). 



/(■^n+1/2) « 
y"('^n+l/2) 



lAh 
1 

2^ 
1 



(y«+2-3j„+i +3j„-);„_i), 



(33) 



where x„+i/2 = ^""^2"^' the scheme's center. 



We will be interested in the behaviour of solutions near singularities (blow up in the first 
derivative). 
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4.1. 5/3(2, R) 

Figure 1 shows the behaviour of the standard and symmetric methods for (20). The choice for 
a corresponds to some since a is the integration constant. The exact solution is ([2T]) and is 
represented by the continuous line. The Newton method applied to the standard scheme con- 
verges poorly (and even fails to converge for h ~ 0.01 and smaller). It stops to describe the 
behaviour of the solution correctly near the first derivative blow up. The symmetric method 
integrates around the circle without any difficulties. 



The mesh equation for the symmetric method is chosen to be 

£2(1",, I'l^'ji^') = -11 = 0. (34) 

This mesh equation is also chosen for the 3rd order equation and for both equations invariant 
under sU. 

Figure 2 shows the behaviour of each method for the 3rd order equation 

x^l + y'^)y"' - x^{y - l)/'2 - lxy'(\ + y'^)y" + y'\\ + y'^f = 0. (35) 

Matlab solver (ode45) encounters a singuarity near x = 1.28. The solution stays finite, 
but, again, there is a blow up in it's first derivative. While ode45 stops integrating and the 
standard method blows up, the symmetric method integrates through the singularity and stays 
finite. 



While the first example is didactic since the analytic solution is known, the second one 
shows that the nice behaviour of the symmetric method holds for the more complex equation 

m. 



Moreover, note the relative simplicity of the symmetric scheme for the 3rd order equation 

(Xn+i(2 + I]) - Xn(2 + /3l))x„+2 + 2(j„+i - j„)j„+2 = 4+1 ' 4 + ' yl (36) 

{X„^2 - (1 + /?/2)^„+l)' + (yn+2 - yn+lf = d + lV'^)^^l4+l 

where jS„ and /i are some constants at each step. The only unknowns in this system are Xn+i 
andjn+2. Thus, solving the system ( [361 ) amounts to finding the intersection between a straight 
line and a circle at each step while the standard finite diff'erence scheme involves a Newton 
iteration on a line counting several hundred characters. Matlab solver, while very precise, also 
has a high computationnal cost. There are known numerical algorithms to find conic intersec- 
tions (we used filOi ). The geometrical similarity between the exact solution for the 2nd order 
equation and this scheme is also interesting. 



4.2. 5/4(2, R) 



Figure 3 shows the behaviour of the standard and symmetric methods for (27). The stan- 
dard method stops correctly describing the solution near the blow up in the first derivative 
(it rapidly diverges after the strange behaviour shown in the figure). The symmetric method 
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9.5 




6.5 1 ' ' ' ' 1 

1 1.5 2 2.5 3 3.5 

X 

Figure 1. Symmetric and standard method for pO] ) with initial conditions {xo - l,yo = 2i,C - 
2,a = 1) 



integrates on the entire branch of the hyperbola without any difficulties. 

Figure 4 shows the behaviour of each method for the 3rd order equation 

2x^(y'^ - l)y"' + (y'^ - l)\8y'^ - 3) + I0xy'y"(y'^ - 1) - xYHSy' - 5) (37) 
= 

Again, the symmetric method integrates through the blow up in the first derivative while 
the other methods stop integrating. 

Similarly as for the 5/3 realization, the standard method for the 3rd order equation leads 
to a complicated nonlinear equation while the symmetric scheme is given by 

{-2(Xn+l — Xn) + A{Xn+iP — Xnqn))Xn+2 + 2(jn+l ~ J«)jn+2 

= yli-yl-(^ly-4) (38) 

where ^„ and p are constants at each step. Solving the symmetric scheme then amounts to 
finding the intersection between a straight line and a hyperbola at each step. There is again 
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Figure 2. All methods for (35 i with initial conditions {xo - l,yo = ^^y'o - IjJo = 3} 



a geometrical similarity between the discrete scheme and the exact solution for the 2nd order 
ODE. 



5. Conclusions 



From this article and two previous ones [jVl [8l we conclude the following. For all 4 inequiv- 
alent realizations of sl{2, R) symmetry preserving discretization provide quahtatively better 
numerical solutions than common finite diff"erence methods (including Matlab's solvers), par- 
ticularly for solutions with singularities. The symmetry preserving methods also provide 
solutions at a lower computational cost. 

We would like to emphasize that we are not simply presenting examples when the 
symmetry preserving numerical methods work. All second order ODEs allowing a nontrivial 
symmetry algebra were discretized in this manner in [l5l|. In fSl] it was shown that if the 
symmetries are Lagrangian ones, then the obtained difference systems can be solved exactly 
(analytically). Article \Jj was devoted to second and third order ODEs with three dimensional 
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Lie algebras (solvable, or simple) and mainly discussed the improved precision of symmetry 
respecting discretizations (at no additional computational cost). The emphasis in [8J and the 
present article is on the behaviour of solutions at singularities (for all 4 distinct realizations of 
sl{2,W)). 
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